#R version 4.4.3 (2025-02-28 ucrt)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 11 x64 

library(dplyr) #1.1.4
library(ggplot2) #3.5.1
library(readr) #2.1.5
library(haven) #2.5.4
library(estimatr) #1.0.6
library(broom.mixed) #0.2.9.6
library(lme4) #1.1-37
library(countrycode) #1.6.1
library(forcats) #1.0.0
library(modelsummary) #2.3.0
library(vdemdata) #15.0
library(clubSandwich) #0.6.0

`%notin%` <- function(x,y) !(x %in% y) 

setwd("C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/cps/replication")

#load data and set up variables for analysis####
gwf <- read_csv("data/gwf_panel.csv") %>% 
  mutate(COWcode = countrycode(countryname, origin = "country.name", destination = "cown")) |> 
  select(countryname, year, COWcode, v2xps_party, v2x_polyarchy, gwf_duration, gwf_regime, gwf_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm) %>% 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.nan(.) | gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.na(.), 0, .)) %>% 
  filter(gwf_prior %in% c("Party", "Military", "Personal") & gwf_regime == "democracy") %>% 
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) %>% 
  mutate(gwf_prior = factor(gwf_prior, levels = c("Party", "Personal", "Military"))) %>% 
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) |> 
  mutate(reactive = -1*ruling)

#read in V-Dem data
vdem <- vdem |> 
  select(country_name, year, COWcode, v2xel_frefair)

#join with gwf data
df <- left_join(gwf, vdem, by = c("COWcode", "year")) |> 
  rename(ccode = COWcode)

gwf1 <- df |> select(ccode, year, v2xel_frefair, v2xps_party, gwf_prior, pi1, prior_pi, reactive, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()

#write in Stata format in order to get PCSEs
write_dta(gwf1, "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/cps/replication/data/stata/stata_replication_table3.dta")

#Models 2 and 4 get PCSE from Stata .do files and replace SEs in text. ##
m1 <- lmer(v2xps_party ~ gwf_prior + pi1 + prior_pi + reactive + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf1, REML=FALSE)
m2 <- lm(v2xps_party ~ gwf_prior + pi1 + prior_pi + reactive + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf1)
m3 <- lmer(v2xel_frefair ~ gwf_prior + pi1 + prior_pi + reactive + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf1, REML=FALSE)
m4 <- lm(v2xel_frefair ~ gwf_prior + pi1 + prior_pi + reactive + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf1)

se_m1 <- vcovCR(m1, cluster = gwf1$ccode, type = "CR2")
se_m2 <- vcovCR(m2, cluster = gwf1$ccode, type = "CR2")
se_m3 <- vcovCR(m3, cluster = gwf1$ccode, type = "CR2")
se_m4 <- vcovCR(m4, cluster = gwf1$ccode, type = "CR2")

sd(residuals(m2))
sd(residuals(m4))

modellist <- list("ML - PI" <- m1, #just baseline
                  "OLS - PI" <- m2,
                  "ML - Elections" <- m3, #just baseline
                  "OLS - Elections" <- m4) #shows basline with time-variant indicators)

cluster_selist <- c(se_m1, se_m2, se_m3, se_m4)

modelsummary(modellist, stars = TRUE, gof_omit = "R2|AIC|RM|F|Log.Lik", robust = cluster_selist, coef_map = c("pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "reactive" = "Reactive ASP",
                                                                                                              "gwf_priorPersonal" = "Personalist", "gwf_priorMilitary" = "Military",
                                                                                                              "gwf_pduration" = "Prior Duration", "regime_elecs" = "Elections", 
                                                                                                              "SD (Intercept ccode)" = "SD: Country", "SD (Observations)" = "SD: Residuals"), output = "latex")